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A boxy bulge in the Milky Way. 
Inversion of the stellar statistics equation with 2MASS data 
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Abstract. Inverting the stellar statistics equation from 2MASS star counts, we obtain the 3D density distribution 
of the Galactic bulge as well as its luminosity function in the K-band. This results in a boxy bulge with axial 
ratios 1:0.5:0.4 and a major axis angle with respect to the Sun-galactic center of 20° — 35°. 

Key words. Galaxy: structure — Infrared: stars 



> ■ 

oo : 

o . 
\o . 

o : 

in . 
o ■ 

^: 

Oh; 
O ' 

a: 



1. Introduction 

Due to extinction, near-infrared (NIR) surveys are the 
most appropriate observations to study the stellar popula- 
tion of the Galactic central components. However, we see 
them edge-on, so the task of constraining their parameters 
is not as easy as in other external galaxies, which can be 
observed face-on. In order to recover the 3D information, 
Lopez-Corredoira et al. (2000; hereafter LOO) developed a 
method of inversion of the star counts based on Bayesian 
statistics, which permits the obtaining of both the stel- 
lar density distribution and the luminosity function (LF) 
without making any a priori assumptions about them. 

The LOO inversion technique has been applied to the 
K-band TMGS star counts up to ttik = 9-0 in some 
strips that crossed the Galactic plane, mostly in positive 
longitudes (LOO). The inversion was also applied to mid- 
infrared MSX Galactic plane maps (Lopez-Corredoira et 
al. 2001a), although in this last case there were very few 
stars to derive the stellar density and only the LF could 
be derived with some accuracy. Now, we have a deeper 
and larger coverage (the whole sky) survey: 2MASS, with 
3 NIR filters (which permits extinction correction) as well 
as some improvements on the disc models (necessary to be 
subtracted from the total star counts and isolate the bulge 
counts), so we have the chance to obtain more precise re- 
sults from the bulge density distribution and its LF. This 
is the goal of the present paper: using the LOO inversion 
technique with the newest data and information. 



2. 2MASS star counts and disc subtraction 

The 2MASS survey provides magnitudes of all stars in 
the NIR filters H, K s in the range 7.5 < tuh < 15.1, 
7.0 < tuk < 14.3 and nearly one magnitude less in the 
galactic plane (Skrutskie et al. 1997). It also provides J 
magnitudes, but we will not use these here. We can build 
maps of cumulative star counts within these ranges (Alard 
2001). In order to make an approximate correction of the 
extinction (Alard 2001), we calculate 



m e =m K - 1.77(H-K), 



(1) 



and we produce maps of cumulative star counts with 
m e < m 6) o and m e ^ between 7.0 and 10.75 in steps of 
0.25 mag. The angular space resolution both in longi- 
tude and latitude is 0.5° in the range -20° < I < 20°, 
-12° < b < 12°. 

In order to isolate the bulge star counts, we must first 
subtract the disc star counts in the regions where both disc 
and bulge contribute along the line of sight. We cannot 
apply the inversion of LOO to disc+bulge counts since both 
components have different LFs and, moreover, the LOO 
method does not work for stars near the Sun as the relative 
error of the kernel in the inversion becomes very large. 

The disc counts are given in Lopez-Corredoira et al. 
(2004), interpreted as a truncated and/or flared disc: 

Disc model 1: An exponential function of the galactocen- 
tric distance (R) for R > 4 kpc: 



( R ~ R Q \ 
\ 1970 pc J , 



1*1 



( R ) star pc 
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h z (R) = 285e 



(lT2 



R—Rq \ 
-0.6K(kpc)] kpc, 



pc, 



(2) 



(3) 
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and a constant density for R < 4 kpc: 

1*1 _o 

p = 5.1p e h ^ R ^ star pc , 
h z (R) = 509 - 48i?(kpc) pc. 



(4) 
(5) 



Disc model 2: A hole in the inner disc which can be mod- 
eled \/R < Rq as: 



p Q e 



(t 



, 3740 pc > 
970 P c"r Rq j 



( R , 3740 pc \ \z\ 
g V 1970 p C T i? ,/g ^-z(H) 



star pc , 

h z (R) = 285[1 + 0.21 kpc _1 (# - R Q ) 
+ 0.056 kpc _2 (#-i? ) 2 ] pc. 



(6) 



In all cases we adopt Rq = 7.9 kpc. Indeed both models 
are almost equivalent for 2.5 kpc< R < Rq. They differ 
substantially only for R < 2.5 kpc. We do not know the 
exact disc profile in this region where it is mixed with the 
bulge, although we know that it lies between model 1 and 
2. 

The disc LF used in the K-band is that given by 
Eaton et al. (1984). In order to obtain the correspon- 
dence between mx and m e , through eq. QJ, we ca ^ _ 
brate ((H — K)q) (color in regions without extinction) 
in a region towards the center where the disc is almost 
isolated and without extinction (I = 20°, b = 6°), as a 
function of m e ,o- In the regions —20.25° < I < —17.75°, 
3.75° < \b\ < 12.25°, we calibrate the amplitude of the 
density which gives the best fit: p© = 0.050 star pc -3 
for both models. This is close to the value used in Lopez- 
Corredoira et al. (2004). 

The results of the disc subtraction using model 1 are 
shown in Fig. E] (Total - disc) for m e < 9.5; for different 
m e ,o the results are similar. 

3. Inversion of bulge star counts 

The fitting of parameters of the density and the LF in 
the bulge would be a way to carry out the analysis, but 
many free parameters would be needed, especially for the 
LF, and a knowledge of the shape of both functions is 
required. The LOO inversion technique is not limited in 
this way. However, some assumptions are needed: the disc 
models in ^J21and the reddening law of eq. (QJ; the assump- 
tion of constant LF throughout the bulge, and that the 
noise effects in the inversion do not significantly change 
the solution (as tested). 

The procedure to invert the stellar statistics equation 
once the disc is subtracted is described in LOO. It is an 
iteration of consecutive Lucy iterations. The LFs in Ri- 
band and the density for each line of sight are obtained as 
numerical functions; nevertheless, we fit "a posteriori" the 
parameters of an analytical 3-D function for the density 
of a triaxial bulge, and we use a weighted average of LFs 
in the different lines of sight. The result of each individual 




Total - disc - bulge 



Fig. 1. Star counts maps with m e < 9.5, for the region 
-20° < I < 20°, -12° < b < 12°. "Total" map stands for 
the counts derived from 2MASS. "Disc" map is the pre- 
diction of the disc "model 1" contribution (see ^2j. "Total- 
disc" is the subtraction of the disc in the total counts map. 
"Bulge" stands for the best triaxial boxy-isodensity fit of 
the result of the inversion (see of "Total-disc" . "Total- 
disc-bulge" are the residuals after subtracting the bulge 
from the "Total-disc" map. The isocounts contours repre- 
sent, from the outermost to the innermost: 600 to 30000 
star deg -2 with steps of 600 star deg -2 . 



iteration does not depend significantly on the initial guess 
and, in the final global outcome, the average/fit of all the 
lines of sights, the dependence on the initial guesses is 
even smaller, much smaller than other sources of error 
and dispersion with respect to the average/fit. 

There is a wrong statement in LOO. Rq cannot be de- 
termined by this method, because any factor in the dis- 
tance scale can be compensated for by making the abso- 
lute magnitudes of all stars fainter or brighter; we recently 
realized that the value Rq = 7.9 kpc given by LOO was in- 
deed obtained because of the selected range of Mk in the 
algorithm. This value was obtained by chance. However, 
the rest of the parameters are independent of the assumed 
value of Rq ; it only affects the horizontal calibration of the 
LF, and, since Rq = 7.9 kpc is in agreement with other 
measures, the presented results in LOO should be correct. 
Thus, we continue to adopt Rq = 7.9 kpc. 
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We repeat here the steps explained by LOO but with 
the present 2MASS data instead of TMGS data. The se- 
lected region is -20° < I < 20°, 2.25° < \b\ < 12.25°; 
we avoid the Galactic plane, where other possible sources 
apart from the disc and bulge may be present. The sub- 
tracted disc model is different from that of LOO (see . 
A major improvement is that we do not need to spec- 
ify the extinction, but we can correct the star counts by 
means of eq. JIJ. By doing this, we avoid the possible 
problems with the patchiness of the extinction. The se- 
lected range of magnitudes to obtain the 3D-density is: 
7.00 < m e ,o < 9.50. For the LF, the inversion is less sen- 
sitive to the errors, so we have further extended the range 
up to m e? o < 10.75. One cannot extend the fainter limit 
too much because the ratio bulge/disc becomes poorer. 
The normalization of the LF is made assuming the value 
given by LOO for <j)(M K = -6.4) = 2.34 x 10~ 5 . We get 
two outcomes: the LF in the range between Mk = — 6.6 
and Mk = —2.8, and the 3D-density distribution of the 
bulge. 

Fig. [21 shows the LF, with the extension to brighter 
stars derived from the TMGS by LOO. In total, we have a 
LF over 8 magnitudes. A comparison with Wainscoat et 
al.'s (1992) synthetic LF (derived from a model of stel- 
lar populations) shows that there is a strong deficit of 
very bright stars {Mk < —8) — supergiants and late M 
giants — with respect to the synthetic LF (LOO), and slight 
differences around Mk = — 6 and Mk = —3.5. 

Fig. shows different cuts of the 3D density distri- 
bution with z =constant, obtained when using model 1 
of the disc. Model 2 of the disc gives similar results (see 
Fig. 0}. This figure gives the smoothed raw result from 
the LOO inversion technique, p(r), and clearly shows the 
non-circularity of the contours in the bulge cuts. This con- 
firms that the bulge is non-axisymmetric, with its elon- 
gated structure closer to us in the first quadrant (I > 0). 

If we assume that the bulge is constituted by triax- 
ial ellipsoids, their parameters would be as follows: For 
the model 1 of the disc, the fitted mean axial ratios are 
1:0.49:0.37; the angle between the major axis of the bulge 
(in the first quadrant) and the line Sun-Galactic center: 
a = 29°; and the density profile p{t) « 6.6e~ 740 p° star 
pc~ 3 , with t = y/x\ + (x 2 /0.49) 2 + (xs/0.37) 2 where x u 
i = 1,2,3 are along the major, middle and minor axis 
of the triaxial ellipsoid, respectively. Tests to recover a 
given 3D-density/LF with known parameters applying the 
whole LOO-method to its projected star counts show that 
typical systematic deviations due only to the inversion 
method and noise are: 90 pc for the scale length of the 
profile, 8° for the angle, and 0.04 for both axial ratios. 
For the model 2 of the disc, the mean axial ratios are 
1:0.49:0.40, angle a = 27°. 

The bulge is closer to axisymmetry (ratios l:l:x) than 
1:0.33:0.22 (Dwek et al. 1995), 1:0.38:0.26 (Freudenreich 
1998) or l:(0.3-0.4):0.3 (Bissantz & Gerhard 2002). 
However, it is not very far from the values of 1:0.6:0.4 
(Binney et al. 1997), 1:0.43:0.29 (Stanek et al. 1997) or 
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Fig. 2. Luminosity function derived in this paper (with 
disc models 1 and 2; error bars correspond to model 1 and 
are similar for model 2) and LOO in comparison with the 
synthetic LF of Wainscoat et al. (1992). Assumed Rq = 
7.9 kpc. 

1:0.54:0.33 (LOO). A lower ellipticity is more in accord with 
other galaxies, because bulges generally do not have large 
asymmetries. This also contributes to clarify the difference 
between a bulge and a bar in our Galaxy. The present 
structure looks like a bulge rather than a stick-like bar. 
We suspect that the derived bulge axial ratios are differ- 
ent due to the dependence of the methods to correct for 
extinction and the disc subtraction. In our case we used 
NIR star counts corrected for extinction based on redden- 
ing, and an inner truncation and/or flare in the disc model 
increases the number of stars associated with the bulge, 
which reduces its global axisymmetry. The result might 
also depend on the method of inversion/fitting and the a 
priori assumptions. 

A parameter which seems to be in agreement with most 
other studies is the orientation of the major axis of the 
bulge (when they refer to the bulge and not to the hypo- 
thetical in-plane bar; see Lopez-Corredoira et al. 2001b): 
12°-30° (Dwek et al. 1995; Nikolaev & Weinberg 1997; 
Stanek et al. 1997; Binney et al. 1997; Freudenreich 1998; 
Sevenster et al. 1999; LOO; Bissantz & Gerhard 2002). 

4. Boxy Bulge 

Although a triaxial ellipsoid represents the bulge in first 
order, the shape of the projection ("Total-disc" in Fig.GJ 
is boxier than the projection of ellipsoids (Fig. [SJ). The 
boxy shape of the Bulge was noted previously (Kent et 
al. 1991, 1992; Dwek et al. 1995). Here, the boxiness is 
even more remarkable than in previous papers, due per- 
haps to a different subtraction of the disc and the avail- 
ability of star counts instead of flux maps. The isoden- 
sity contours in Fig. 03 especially those at higher t and 
smaller z, look boxier also than an ellipse. If we use 
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Fig. 3. Different cuts (z = ±440, ±770, ±1100 pc) 
through the 3D bulge density distribution, obtained from 
the bulge counts with the LOO inversion technique af- 
ter subtracting the model 1 for the disc. Vertical axis, 
parallel to the line Sun-Galactic center, between 4.4 and 
11.0 kpc from the Sun; increased distance from the Sun 
upwards. Horizontal axis, perpendicular to the line Sun- 
Galactic center: between -2.2 and 2.2 kpc; positive lon- 
gitude to the left. Rq = 7.9 kpc. The contours repre- 
sent isodensity regions; the density range between and 
1.4 star pc -3 . Gaussian smoothing of the data was ap- 
plied with a x = (j y = 220 pc, a 2 (±440 pc) = 55 pc, 
^(±770 pc) = 110 pc, cr 2 (±1100 pc) = 110 pc. 




Fig. 4. Left: Star counts of the disc "model 2" contribu- 
tion (see ®; isocounts contours like in Fig. E] Middle: 
"Total-disc" with "model 2". Right: Cut at z = -440 pc 
as in Fig. 01 but with "model 2" for the disc. 



t = [x A ± (x2/a) 4 ± (X3/6) 4 ] 1 / 4 — a similar boxy structure 
as was proposed by Kent et al. (1991) for the axisymmet- 
ric case — we get a better representation with somewhat 
fewer residuals (Fig. 0. For this boxy- model, we adopt 
the same angle a and axial ratios, and the same profile 
p(t) as for the ellipsoids but with tboxy = 0.866t e iiip- This 
is the approximation that gives the best correspondence 
between ellipsoids and boxy-isocontours with the same in- 
clination and axial ratios. The integrated boxy-contours 
do not perfectly represent the observed isocount contours 
in the outmost regions — it seems that the corners are not 
so abrupt — but they give a better model than the ellipti- 
cal bulge. The asymmetry of the counts in the residuals 



in plane regions (Fig. d "total-disc-bulge" ) suggests that 
part of these remaining stars are due to another axisym- 
metric component, presumably a long bar. Other compo- 
nents like a ring, spiral arms, and errors in the disc sub- 
traction and extinction correction in plane regions are also 
expected to contribute to the residuals. 




Fig. 5. Left: best triaxial ellipsoid fit of the result of the 
inversion (see © of "Total-disc" in Fig.ffl Right: "Total- 
disc-bulge" map. Isocount contours as in Fig. [I] 



5. Discussion and conclusions 

It is well known that boxy bulges are related to bars 
in edge-on galaxies (Merrifield & Kuijken 1999). Their 
bar-like nature is revealed by the spectroscopical data 
which show the typical dynamics of multiple-component 
line emission with one of the components being bar-like. 
Our Galaxy has one of these boxy-bulges, which is not 
itself long-bar-like but quite thick: axial ratios 1:0.5:0.4. 
Therefore, the scenario of an in-plane bar apart from the 
bulge (Lopez-Corredoira et al. 2001b, Hammersley et al. 
2001, and references therein) is supported. The asymme- 
try of the residuals "Total-disc-bulge" also agrees with this 
scenario. The triaxial bulge has an angle a « 20° — 35° 
with respect to the line Galactic center-Sun, smaller than 
the long-bar angle. 
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